#include "LineAbDipole.h"
#include "mfAddFuncCalcTask.h"
#include "BesselFunctions.h"
#include "kernel/feObserver.h"
#include "kernel/feDefs.h"
#include "FastMath.h"
#include <math.h>
#include "kernel/fePolylineGenerator.h"
#include "kernel/feMaterial.h"
#include "Timer.h"
#include <omp.h>

//---------------------------------------------------------------------------------------------------------------------------------------------
const InfoM  LineAbDipole::InfosM[194] = {
	{      1.000000e-007  ,  1.000e-003  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.154780e-007  ,  1.000e-003  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.333520e-007  ,  1.000e-003  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.539930e-007  ,  1.000e-003  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.778280e-007  ,  9.000e-004  ,  9.000e+007  ,  1.050e+000  } ,
	{      2.053530e-007  ,  8.000e-004  ,  8.000e+007  ,  1.050e+000  } ,
	{      2.371370e-007  ,  7.000e-004  ,  7.000e+007  ,  1.050e+000  } ,
	{      2.738420e-007  ,  6.000e-004  ,  6.000e+007  ,  1.050e+000  } ,
	{      3.162280e-007  ,  5.000e-004  ,  5.000e+007  ,  1.050e+000  } ,
	{      3.651740e-007  ,  5.000e-004  ,  5.000e+007  ,  1.050e+000  } ,
	{      4.216970e-007  ,  4.000e-004  ,  4.000e+007  ,  1.050e+000  } ,
	{      4.869680e-007  ,  3.000e-004  ,  3.000e+007  ,  1.050e+000  } ,
	{      5.623410e-007  ,  2.000e-004  ,  2.000e+007  ,  1.050e+000  } ,
	{      6.493820e-007  ,  1.000e-004  ,  1.000e+007  ,  1.050e+000  } ,
	{      7.498940e-007  ,  1.000e-004  ,  1.000e+007  ,  1.050e+000  } ,
	{      8.659640e-007  ,  1.000e-004  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.000000e-006  ,  1.000e-004  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.154780e-006  ,  1.000e-004  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.333520e-006  ,  1.000e-004  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.539930e-006  ,  1.000e-004  ,  1.000e+007  ,  1.050e+000  } ,
	{      1.778280e-006  ,  9.000e-005  ,  9.000e+006  ,  1.050e+000  } ,
	{      2.053530e-006  ,  8.000e-005  ,  8.000e+006  ,  1.050e+000  } ,
	{      2.371370e-006  ,  7.000e-005  ,  7.000e+006  ,  1.050e+000  } ,
	{      2.738420e-006  ,  6.000e-005  ,  6.000e+006  ,  1.050e+000  } ,
	{      3.162280e-006  ,  5.000e-005  ,  5.000e+006  ,  1.050e+000  } ,
	{      3.651740e-006  ,  5.000e-005  ,  5.000e+006  ,  1.050e+000  } ,
	{      4.216970e-006  ,  4.000e-005  ,  4.000e+006  ,  1.050e+000  } ,
	{      4.869680e-006  ,  3.000e-005  ,  3.000e+006  ,  1.050e+000  } ,
	{      5.623410e-006  ,  2.000e-005  ,  2.000e+006  ,  1.050e+000  } ,
	{      6.493820e-006  ,  1.000e-005  ,  1.000e+006  ,  1.050e+000  } ,
	{      7.498940e-006  ,  1.000e-005  ,  1.000e+006  ,  1.050e+000  } ,
	{      8.659640e-006  ,  1.000e-005  ,  1.000e+006  ,  1.050e+000  } ,
	{      1.000000e-005  ,  1.000e-005  ,  1.000e+006  ,  1.050e+000  } ,
	{      1.154780e-005  ,  1.000e-005  ,  1.000e+006  ,  1.050e+000  } ,
	{      1.333520e-005  ,  1.000e-005  ,  1.000e+006  ,  1.050e+000  } ,
	{      1.539930e-005  ,  1.000e-005  ,  1.000e+006  ,  1.050e+000  } ,
	{      1.778280e-005  ,  9.000e-006  ,  9.000e+005  ,  1.050e+000  } ,
	{      2.053530e-005  ,  8.000e-006  ,  8.000e+005  ,  1.050e+000  } ,
	{      2.371370e-005  ,  7.000e-006  ,  7.000e+005  ,  1.050e+000  } ,
	{      2.738420e-005  ,  6.000e-006  ,  6.000e+005  ,  1.050e+000  } ,
	{      3.162280e-005  ,  5.000e-006  ,  5.000e+005  ,  1.050e+000  } ,
	{      3.651740e-005  ,  5.000e-006  ,  5.000e+005  ,  1.050e+000  } ,
	{      4.216970e-005  ,  4.000e-006  ,  4.000e+005  ,  1.050e+000  } ,
	{      4.869680e-005  ,  3.000e-006  ,  3.000e+005  ,  1.050e+000  } ,
	{      5.623410e-005  ,  2.000e-006  ,  2.000e+005  ,  1.050e+000  } ,
	{      6.493820e-005  ,  1.000e-006  ,  1.000e+005  ,  1.050e+000  } ,
	{      7.498940e-005  ,  1.000e-006  ,  9.000e+004  ,  1.050e+000  } ,
	{      8.659640e-005  ,  1.000e-006  ,  8.000e+004  ,  1.050e+000  } ,
	{      1.000000e-004  ,  1.000e-006  ,  7.000e+004  ,  1.050e+000  } ,
	{      1.154780e-004  ,  9.000e-007  ,  6.000e+004  ,  1.050e+000  } ,
	{      1.333520e-004  ,  8.000e-007  ,  5.000e+004  ,  1.050e+000  } ,
	{      1.539930e-004  ,  7.000e-007  ,  4.000e+004  ,  1.050e+000  } ,
	{      1.778280e-004  ,  5.000e-007  ,  3.000e+004  ,  1.050e+000  } ,
	{      2.053530e-004  ,  4.000e-008  ,  2.000e+004  ,  1.050e+000  } ,
	{      2.371370e-004  ,  3.000e-008  ,  1.000e+004  ,  1.050e+000  } ,
	{      2.738420e-004  ,  2.000e-008  ,  9.000e+003  ,  1.050e+000  } ,
	{      3.162280e-004  ,  1.000e-008  ,  8.000e+003  ,  1.050e+000  } ,
	{      3.651740e-004  ,  1.000e-009  ,  7.000e+003  ,  1.050e+000  } ,
	{      4.216970e-004  ,  1.000e-009  ,  6.000e+003  ,  1.050e+000  } ,
	{      4.869680e-004  ,  1.000e-009  ,  5.000e+003  ,  1.050e+000  } ,
	{      5.623410e-004  ,  1.000e-009  ,  4.000e+003  ,  1.050e+000  } ,
	{      6.493820e-004  ,  1.000e-010  ,  3.000e+003  ,  1.050e+000  } ,
	{      7.498940e-004  ,  1.000e-010  ,  2.000e+003  ,  1.050e+000  } ,
	{      8.659640e-004  ,  1.000e-010  ,  1.000e+003  ,  1.050e+000  } ,
	{      1.000000e-003  ,  1.000e-011  ,  8.000e+002  ,  1.050e+000  } ,
	{      1.154782e-003  ,  1.000e-011  ,  8.000e+002  ,  1.050e+000  } ,
	{      1.333521e-003  ,  1.000e-011  ,  7.000e+002  ,  1.050e+000  } ,
	{      1.539927e-003  ,  1.000e-011  ,  6.000e+002  ,  1.050e+000  } ,
	{      1.778279e-003  ,  1.000e-011  ,  6.000e+002  ,  1.050e+000  } ,
	{      2.053525e-003  ,  1.000e-011  ,  5.000e+002  ,  1.050e+000  } ,
	{      2.371374e-003  ,  1.000e-011  ,  4.000e+002  ,  1.050e+000  } ,
	{      2.738420e-003  ,  1.000e-011  ,  4.000e+002  ,  1.050e+000  } ,
	{      3.162278e-003  ,  1.000e-011  ,  3.000e+002  ,  1.050e+000  } ,
	{      3.651741e-003  ,  1.000e-011  ,  2.000e+002  ,  1.050e+000  } ,
	{      4.216965e-003  ,  1.000e-011  ,  2.000e+002  ,  1.050e+000  } ,
	{      4.869675e-003  ,  1.000e-011  ,  1.000e+002  ,  1.050e+000  } ,
	{      5.623413e-003  ,  1.000e-011  ,  9.000e+001  ,  1.050e+000  } ,
	{      6.493816e-003  ,  1.000e-011  ,  8.000e+001  ,  1.050e+000  } ,
	{      7.498942e-003  ,  1.000e-011  ,  7.000e+001  ,  1.050e+000  } ,
	{      8.659643e-003  ,  1.000e-012  ,  6.000e+001  ,  1.050e+000  } ,
	{      1.000000e-002  ,  1.000e-012  ,  5.000e+001  ,  1.050e+000  } ,
	{      1.154782e-002  ,  1.000e-012  ,  4.000e+001  ,  1.050e+000  } ,
	{      1.333521e-002  ,  1.000e-012  ,  3.000e+001  ,  1.050e+000  } ,
	{      1.539927e-002  ,  1.000e-012  ,  2.000e+001  ,  1.050e+000  } ,
	{      1.778279e-002  ,  1.000e-012  ,  1.000e+001  ,  1.050e+000  } ,
	{      2.053525e-002  ,  9.000e-012  ,  1.000e+001  ,  1.050e+000  } ,
	{      2.371374e-002  ,  8.000e-012  ,  1.000e+001  ,  1.050e+000  } ,
	{      2.738420e-002  ,  7.000e-012  ,  9.000e+000  ,  1.050e+000  } ,
	{      3.162278e-002  ,  6.000e-012  ,  8.000e+000  ,  1.050e+000  } ,
	{      3.651741e-002  ,  5.000e-012  ,  8.000e+000  ,  1.050e+000  } ,
	{      4.216965e-002  ,  4.000e-012  ,  7.000e+000  ,  1.050e+000  } ,
	{      4.869675e-002  ,  3.000e-012  ,  6.000e+000  ,  1.050e+000  } ,
	{      5.623413e-002  ,  3.000e-012  ,  6.000e+000  ,  1.050e+000  } ,
	{      6.493816e-002  ,  2.000e-012  ,  5.000e+000  ,  1.050e+000  } ,
	{      7.498942e-002  ,  2.000e-012  ,  5.000e+000  ,  1.050e+000  } ,
	{      8.659643e-002  ,  1.000e-012  ,  5.000e+000  ,  1.050e+000  } ,
	{      1.000000e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{      1.074608e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{      1.154782e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{      1.240938e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.333521e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.433013e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.539927e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.654817e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.778279e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.910953e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     2.053525e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     2.206734e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     2.371374e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     2.548297e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     2.738420e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     2.942727e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     3.162278e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     3.398208e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     3.651741e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     3.924190e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     4.216965e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     4.531584e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     4.869675e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     5.232991e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     5.623413e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     6.042964e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     6.493816e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     6.978306e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     7.498942e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     8.058422e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     8.659643e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     9.305720e-001  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.000000e+000  ,  7.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.074608e+000  ,  6.000e-013  ,  1.000e+000  ,  1.050e+000  } ,
	{     1.154782e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     1.240938e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     1.333521e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     1.433013e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     1.539927e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     1.654817e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     1.778279e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     1.910953e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     2.053525e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     2.206734e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     2.371374e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     2.548297e+000  ,  6.000e-013  ,  5.000e-001  ,  1.050e+000  } ,
	{     2.738420e+000  ,  6.000e-013  ,  1.000e-001  ,  1.050e+000  } ,
	{     2.942727e+000  ,  6.000e-013  ,  1.000e-001  ,  1.050e+000  } ,
	{     3.162278e+000  ,  6.000e-013  ,  1.000e-001  ,  1.050e+000  } ,
	{     3.398208e+000  ,  6.000e-013  ,  1.000e-001  ,  1.050e+000  } ,
	{     3.651741e+000  ,  6.000e-013  ,  1.000e-001  ,  1.050e+000  } ,
	{     3.924190e+000  ,  6.000e-013  ,  5.000e-002  ,  1.050e+000  } ,
	{     4.216965e+000  ,  6.000e-013  ,  5.000e-002  ,  1.050e+000  } ,
	{     4.531584e+000  ,  6.000e-013  ,  5.000e-002  ,  1.050e+000  } ,
	{     4.869675e+000  ,  6.000e-013  ,  5.000e-002  ,  1.050e+000  } ,
	{     5.232991e+000  ,  5.000e-013  ,  5.000e-002  ,  1.050e+000  } ,
	{     5.623413e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     6.042964e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     6.493816e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     6.978306e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     7.498942e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     8.058422e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     8.659643e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     9.305720e+000  ,  5.000e-013  ,  1.000e-002  ,  1.050e+000  } ,
	{     1.000000e+001  ,  5.000e-013  ,  1.000e-003  ,  1.050e+000  } ,
	{     1.074608e+001  ,  4.000e-013  ,  1.000e-003  ,  1.050e+000  } ,
	{     1.154782e+001  ,  4.000e-013  ,  1.000e-003  ,  1.050e+000  } ,
	{     1.240938e+001  ,  3.000e-013  ,  1.000e-003  ,  1.050e+000  } ,
	{     1.333521e+001  ,  3.000e-013  ,  1.000e-003  ,  1.050e+000  } ,
	{     1.433013e+001  ,  2.000e-013  ,  1.000e-003  ,  1.050e+000  } ,
	{     1.539927e+001  ,  2.000e-013  ,  1.000e-003  ,  1.050e+000  } ,
	{     1.654817e+001  ,  1.000e-013  ,  5.000e-004  ,  1.050e+000  } ,
	{     1.778279e+001  ,  1.000e-013  ,  5.000e-004  ,  1.050e+000  } ,
	{     1.910953e+001  ,  9.000e-014  ,  5.000e-004  ,  1.050e+000  } ,
	{     2.053525e+001  ,  9.000e-014  ,  5.000e-004  ,  1.050e+000  } ,
	{     2.206734e+001  ,  8.000e-014  ,  5.000e-004  ,  1.050e+000  } ,
	{     2.371374e+001  ,  8.000e-014  ,  5.000e-004  ,  1.050e+000  } ,
	{     2.548297e+001  ,  7.000e-014  ,  1.000e-005  ,  1.050e+000  } ,
	{     2.738420e+001  ,  7.000e-014  ,  1.000e-005  ,  1.050e+000  } ,
	{     2.942727e+001  ,  6.000e-014  ,  1.000e-005  ,  1.050e+000  } ,
	{     3.162278e+001  ,  6.000e-014  ,  1.000e-005  ,  1.050e+000  } ,
	{     3.398208e+001  ,  5.000e-014  ,  1.000e-005  ,  1.050e+000  } ,
	{     3.651741e+001  ,  5.000e-014  ,  1.000e-005  ,  1.050e+000  } ,
	{     3.924190e+001  ,  4.000e-014  ,  5.000e-006  ,  1.050e+000  } ,
	{     4.216965e+001  ,  4.000e-014  ,  5.000e-006  ,  1.050e+000  } ,
	{     4.531584e+001  ,  3.000e-014  ,  5.000e-006  ,  1.050e+000  } ,
	{     4.869675e+001  ,  3.000e-014  ,  5.000e-006  ,  1.050e+000  } ,
	{     5.232991e+001  ,  2.000e-014  ,  5.000e-006  ,  1.050e+000  } ,
	{     5.623413e+001  ,  2.000e-014  ,  5.000e-006  ,  1.050e+000  } ,
	{     6.042964e+001  ,  1.000e-014  ,  5.000e-006  ,  1.050e+000  } ,
	{     6.493816e+001  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  } ,
	{     6.978306e+001  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  } ,
	{     7.498942e+001  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  } ,
	{     8.058422e+001  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  } ,
	{     8.659643e+001  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  } ,
	{     9.305720e+001  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  } ,
	{     1.000000e+002  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  } ,
	{     1.000000e+005  ,  1.000e-014  ,  1.000e-006  ,  1.050e+000  }
};


//---------------------------------------------------------------------------------------------------------------------------------------------
LineAbDipole::LineAbDipole()
{
	MsCount = 0;
	MsArray = 0;

	Layers = 0;
	NumLayers = 0;

	Length = 1;
	SegmentsNum = 0;

	V = 0;
	X = 0;
	DxDt = 0;
	DxDz = 0;
	DvDz = 0;
	NeedRecalc = false;

	TimesCount = 240;
	double tau = 1.1;
	double ftime = 1e-10;
	double step = ftime;
	double t = ftime;	

	//fill time array
	TimesArray = new double[TimesCount];
	for(int i=0; i<TimesCount; i++)
	{
		TimesArray[i] = t;
		step *= tau;
		t += step;
	}
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void LineAbDipole::TimeInterpolation(vector<double> &iFrom, double *oTo, double *iTimes, int iTimeCnt)
{
	//interpolation from main time mesh to extern mesh
	int it = 0;
	for(int i=1; i<TimesCount-1 && it < iTimeCnt;)
	{
		double t = iTimes[it];
		if(TimesArray[i] > t)
		{			
			double ts[3] = {TimesArray[i-1], TimesArray[i], TimesArray[i+1] };
			double fs[3] = { iFrom[i-1], iFrom[i], iFrom[i+1] };
			oTo[it] = feInterpolation::Quadratic(ts, fs, t);
			it++;
		}
		else
			i++;
	}
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double LineAbDipole::J0Int(double iF1, double iF2, double iR, double iX0, double iX1)
{
	BesselFunctions& bf = BesselFunctions::GetInstance();


	double a = iF1/(iX0-iX1)-iF2/(iX0-iX1);
	double b = -iX1*iF1/(iX0-iX1)+iX0*iF2/(iX0-iX1);


	double j0x0 = bf.J0(iX0*iR);
	double j1x0 = bf.J1(iX0*iR);
	double s0x0 = bf.StruveH0(iX0*iR);
	double s1x0 = bf.StruveH1(iX0*iR);

	double j0x1 = bf.J0(iX1*iR);
	double j1x1 = bf.J1(iX1*iR);
	double s0x1 = bf.StruveH0(iX1*iR);
	double s1x1 = bf.StruveH1(iX1*iR);


	double res = -0.5*(2*a*iX0*j1x0 + 2*b*iR*iX0*j0x0 + b*PI*iR*iX0*s0x0*j1x0
		- b*PI*iR*iX0*s1x0*j0x0 - 2*a*iX1*j1x1 - 2*b*iR*iX1*j0x1 -
		b*PI*iR*iX1*s0x1*j1x1 + b*PI*iR*iX1*s1x1*j0x1)/iR;
	return res;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
int LineAbDipole::GetSegmentsNum(double iX, double iY)
{	return 10;
	//get vector coords, y-coords of dipole= 0
	double vx1  = -Length/2 - iX;
	double vy1  = -iY;

	double vx2 = Length/2 - iX;
	double vy2 = -iY;
	
	double cosv = (vx1 * vx2 + vy1 * vy2)/(FastMath::Sqrt((vx1*vx1 + vy1*vy1)*(vx2*vx2 + vy2*vy2)));

	//one segment per 20 degrees
	int segments = (int)(acos(cosv)*180/PI);
	if( segments <= 0 ) segments = 1;
	return segments;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double LineAbDipole::J1Int(double iF1, double iF2, double iR, double iX0, double iX1)
{
	BesselFunctions& bf = BesselFunctions::GetInstance();


	double a = iF1/(iX0-iX1)-iF2/(iX0-iX1);
	double b = -iX1*iF1/(iX0-iX1)+iX0*iF2/(iX0-iX1);

	double j0x0 = bf.J0(iX0*iR);
	double j1x0 = bf.J1(iX0*iR);
	double s0x0 = bf.StruveH0(iX0*iR);
	double s1x0 = bf.StruveH1(iX0*iR);

	double j0x1 = bf.J0(iX1*iR);
	double j1x1 = bf.J1(iX1*iR);
	double s0x1 = bf.StruveH0(iX1*iR);
	double s1x1 = bf.StruveH1(iX1*iR);


	double t1 = a*0.3141592653589793e1;
	double t3 = j1x0;
	double t5 = s0x0;
	double t8 = j0x0;
	double t10 = s1x0;
	
	double t16 = j1x1;
	double t18 = s0x1;
	double t21 = j0x1;
	double t23 = s1x1;

	double res = -(t1*iX0*t3*t5-t1*iX0*t8*t10-2.0*b*t8-t1*iX1*t16*t18+t1*iX1*t21*t23
		+2.0*b*t21)/iR/2.0;
	return res;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double LineAbDipole::ExIntFunction(int iMsIndex, int iTimeIndex, double iR, double iX, double iY, double iSigma)
{
	double m = MsArray[iMsIndex];
	BesselFunctions &bs = BesselFunctions::GetInstance();
	double j0 = bs.J0(m*iR);
	double j1 = bs.J1(m*iR);


	double s0 = -(j0-j1/(m*iR))*m/(iR*iR);
	double s1 = j1/(iR*iR*iR);
	double s2 = -j1/iR;
	return (s0*iX*iX+s1*iX*iX+s2)*DvDz[iTimeIndex][iMsIndex] + (s0*iY*iY+s1*iY*iY+s2)*MU0*DxDt[iTimeIndex][iMsIndex];
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double LineAbDipole::EyIntFunction( int iMsIndex, int iTimeIndex, double iR, double iX, double iY )
{
	double m = MsArray[iMsIndex];
	BesselFunctions &bs = BesselFunctions::GetInstance();
	double j0 = bs.J0(m*iR);
	double j1 = bs.J1(m*iR);


	double s0 = -(j0-j1/(m*iR))*m/(iR*iR);
	double s1 = j1/(iR*iR*iR);

	return (s0+s1)*iX*iY*(DvDz[iTimeIndex][iMsIndex] - MU0 * DxDt[iTimeIndex][iMsIndex]);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double LineAbDipole::EzIntFunction(int iMsIndex, int iTimeIndex, double iR, double iX)
{
	double m = MsArray[iMsIndex];
	BesselFunctions &bs = BesselFunctions::GetInstance();
	double j1 = bs.J1(m*iR);
	//done
	return m*iX/iR*j1*m*V[iTimeIndex][iMsIndex];
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double LineAbDipole::DBzDtIntFunction(int iMsIndex, int iTimeIndex, double iR, double iX, double iY )
{
	return -MsArray[iMsIndex]*MsArray[iMsIndex]*iY*DxDt[iTimeIndex][iMsIndex]/iR;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double LineAbDipole::ByIntFunction(int iMsIndex, double iSigma, int iTimeIndex, double iR, double iX, double iY)
{
	double m = MsArray[iMsIndex];
	BesselFunctions &bs = BesselFunctions::GetInstance();
	double j0 = bs.J0(m*iR);
	double j1 = bs.J1(m*iR);


	double s0 = -(j0-j1/(m*iR))*m/(iR*iR);
	double s1 = j1/(iR*iR*iR);
	double s2 = -j1/iR;
	return (s0*iX*iX+s1*iX*iX+s2)*V[iTimeIndex][iMsIndex]*iSigma + (s0*iY*iY+s1*iY*iY+s2)*DxDz[iTimeIndex][iMsIndex];
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void LineAbDipole::SetMedium(ConductiveLayer *iLayers, int iNumLayers)
{
	NumLayers = iNumLayers;
	Layers = iLayers;
	Zhelper.Load();
	Vhelper.Load();
	NeedRecalc = true;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
bool LineAbDipole::CreateX(double iZ)
{
	if(!Layers)
		return false;

	//delete m
	if(MsArray)
		delete [] MsArray;
	MsCount = 0;

	//m-mesh params
	MsCount = 100;

	double pointsperdegree = 11;
	double startdegree = -6;

	//allocate ms
	MsArray = new double[MsCount];

	//allocate dzdts
	if(!X)
	{
		X = new double *[TimesCount];
		V = new double *[TimesCount];
		DvDz = new double *[TimesCount];
		DxDt = new double *[TimesCount];
		DxDz = new double *[TimesCount];
		for(int it=0; it<TimesCount; it++)
		{
			X[it] = new double[MsCount];
			V[it] = new double[MsCount];
			DvDz[it] = new double[MsCount];
			DxDt[it] = new double[MsCount];
			DxDz[it] = new double[MsCount];
		}
	}
	

#pragma omp parallel shared(pointsperdegree, startdegree) num_threads(2) 
	{	
#pragma omp sections
		{
#pragma omp section
			{
				double m = 1e-6;
				int iinf = 0;

				mfAddFuncCalcTask xtask;
				for(int im=0, i=0; im<MsCount; im++, i++)
				{
					m = pow(10, startdegree+i/pointsperdegree);

					MsArray[im] = m;

					//read task params from table
					while(InfosM[iinf].M < m) iinf++;
					double firstt = max(InfosM[iinf].Tmin, TimesArray[0]);
					double maxt = min(InfosM[iinf].Tmax, TimesArray[TimesCount-1]);
					double timecoef = InfosM[iinf].Tau;

					//get num layers over geometric progression sum
					int nl = (int)(log(1-maxt/firstt * (1-timecoef))/log(timecoef) + 0.5);
					Zhelper.Init(Layers, m);

					//init task
					if(!xtask.Init(Layers, NumLayers, timecoef, firstt, nl, m, &Zhelper,true)) continue;


					//solve
					if(!xtask.Solve()) continue;


					//fill polyline for current m
					//fill z arrays for second task
					for(int it=0; it<TimesCount; it++)
					{
						//X[it][im] = xtask.GetU(iZ, TimesArray[it])/(2*m);
						DxDt[it][im] = xtask.GetDuDt(iZ, TimesArray[it])/(2*m);
						DxDz[it][im] = -xtask.GetDuDz(iZ, TimesArray[it])/(2*m);
					}
				}
			}
#pragma omp section
			{
				double m = 1e-6;
				int iinf = 0;

				mfAddFuncCalcTask vtask;
				for(int im=0, i=0; im<MsCount; im++, i++)
				{
					m = pow(10, startdegree+i/pointsperdegree);

					//read task params from table
					while(InfosM[iinf].M < m) iinf++;
					double firstt = max(InfosM[iinf].Tmin, TimesArray[0]);
					double maxt = min(InfosM[iinf].Tmax, TimesArray[TimesCount-1]);
					double timecoef = InfosM[iinf].Tau;

					//get num layers over geometric progression sum
					int nl = (int)(log(1-maxt/firstt * (1-timecoef))/log(timecoef) + 0.5);
					Vhelper.Init(Layers, m);
					if(!vtask.Init(Layers, NumLayers, timecoef, firstt, nl, m, &Vhelper, true)) continue;


					//solve
					if(!vtask.Solve()) continue;

					for(int it=0; it<TimesCount; it++)
					{
						V[it][im] = vtask.GetU(iZ, TimesArray[it])/Layers[0].Sigma;
						DvDz[it][im] = -vtask.GetDuDz(iZ, TimesArray[it])/Layers[0].Sigma;
					}
				}
			}
		}
	}

	NeedRecalc = false;

// 	fePolylineGenerator g;
// 	double *f = new double[MsCount];
// 	for(int i=0; i<MsCount; i++)
// 	{
// 		f[i] = ExIntFunction(MsArray[i], 120, 5, 3, 4, 1);	
// 	}
// 
// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(MsArray, f, MsCount, 1), "v_0_0001");
// 
// 	for(int i=0; i<MsCount; i++)
// 	{
// 		f[i] = ExIntFunction(MsArray[i], 144, 5, 3, 4, 1);	
// 	}
// 	
// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(MsArray, f, MsCount, 1), "v_0_001");
// 
// 	for(int i=0; i<MsCount; i++)
// 	{
// 		f[i] = ExIntFunction(MsArray[i], 169, 5, 3, 4, 1);	
// 	}
// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(MsArray, f, MsCount, 1), "v_0_01");
// 
// 	
// 	for(int i=0; i<MsCount; i++)
// 	{
// 		f[i] = ExIntFunction(MsArray[i], 193, 5, 3, 4, 1);	
// 	}
// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(MsArray, f, MsCount, 1), "v_0_1");

// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(MsArray, DxDt[100], MsCount, 1), "x_0_0001");
// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(MsArray, DxDt[169], MsCount, 1), "x_0_01");
// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(MsArray, DxDt[193], MsCount, 1), "x_0_1");
	return true;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool LineAbDipole::Ex(double iX, double iY, double iZ, double *oResV, double *iTimes, int iTimeCount)
{
	if(!Layers) return false;
	if(!V || NeedRecalc) CreateX(iZ);

	std::vector<double> ey;
	ey.reserve(TimesCount);

	


	//try to extract from cache
// 	if(EyCache.GetV(iX, iY, iZ, 0.01*r, ey))
// 	{
// 		TimeInterpolation(ey, oResV, iTimes, iTimeCount);
// 		return true;
// 	}

	int n = GetSegmentsNum(iX, iY);
	double step = Length/n;

	for(int it = 0; it<TimesCount; it++)
	{
		double s = 0;
		for(int sn = 0; sn<n; sn++)
		{
			double sx = -Length/2 + step * sn + step/2;
			double x = iX + sx;
			double r = FastMath::Sqrt(x*x + iY*iY);

			//numerical			
			double f[2];
			double m[2];
			f[0] = ExIntFunction(0, it, r, x, iY, Layers[0].Sigma);
			m[0] = MsArray[0];

			for(int i=1; i<MsCount; i++)
			{
				//fill arrays
				f[1] = ExIntFunction(i, it, r, x, iY, Layers[0].Sigma);
				m[1] = MsArray[i];

				s += 0.5*(f[0]+f[1])*(m[1]-m[0]);	

				//update f[0]
				f[0] = f[1];
				m[0] = m[1];			
			}			
		}

		//push solvation
		s *= step;
		s /= 2*PI;
		ey.push_back(s);
	}

/*	EyCache.PushV(iX, iY, iZ, ey);*/
	TimeInterpolation(ey, oResV, iTimes, iTimeCount);

// 	fePolylineGenerator g;
// 	
// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(iTimes, oResV, iTimeCount, 0.1), "c:\\anal");
	return true;	
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool LineAbDipole::Ey( double iX, double iY, double iZ, double *oResV, double *iTimes, int iTimeCount )
{
	if(!Layers) return false;
	if(!V || NeedRecalc) CreateX(iZ);

	std::vector<double> ey;
	ey.reserve(TimesCount);


	//try to extract from cache
	// 	if(EyCache.GetV(iX, iY, iZ, 0.01*r, ey))
	// 	{
	// 		TimeInterpolation(ey, oResV, iTimes, iTimeCount);
	// 		return true;
	// 	}


	int n = GetSegmentsNum(iX, iY);
	double step = Length/n;

	double zero_t = 0.1*Layers[0].H*Layers[0].H  * MU0 * Layers[0].Sigma;
	for(int it = 0; it<TimesCount; it++)
	{
		double s = 0;
		if(TimesArray[it] >= zero_t)
		{
			for(int sn = 0; sn<n; sn++)
			{
				double sx = -Length/2 + step * sn + step/2;
				double x = iX + sx;
				double r = FastMath::Sqrt(x*x + iY*iY);


				//numerical			
				double f[2];
				double m[2];
				f[0] = EyIntFunction(0, it, r, x, iY);
				m[0] = MsArray[0];

				for(int i=1; i<MsCount; i++)
				{
					//fill arrays
					f[1] = EyIntFunction(i, it, r, x, iY);
					m[1] = MsArray[i];

					s += 0.5*(f[0]+f[1])*(m[1]-m[0]);	

					//update f[0]
					f[0] = f[1];
					m[0] = m[1];			
				}			
			}
		}

		//push solvation
		s *= step;
		s /= -2*PI;
		ey.push_back(s);
	}

	/*	EyCache.PushV(iX, iY, iZ, ey);*/
	TimeInterpolation(ey, oResV, iTimes, iTimeCount);

	// 	fePolylineGenerator g;
	// 	
	// 	g.SavePolylineToFile(&g.CreatePolylineFromArray(iTimes, oResV, iTimeCount, 0.1), "c:\\anal");
	return true;	
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool LineAbDipole::Ez(double iX, double iY, double iZ, double *oResV, double *iTimes, int iTimeCount)
{	
	if(!Layers) return false;
	if(!V || NeedRecalc) CreateX(iZ);

	std::vector<double> ez;
	ez.reserve(TimesCount);


	//try to extract from cache
// 	if(EzCache.GetV(iX, iY, iZ, 0.01*r, ez))
// 	{
// 		TimeInterpolation(ez, oResV, iTimes, iTimeCount);
// 		return true;
// 	}

	
	int n = GetSegmentsNum(iX, iY);
	double step = Length/n;

	for(int it = 0; it<TimesCount; it++)
	{
		double s = 0;
		for(int sn = 0; sn<n; sn++)
		{
			double sx = -Length/2 + step * sn + step/2;
			double x = iX + sx;
			double r = FastMath::Sqrt(x*x + iY*iY);


			//numerical			
			double f[2];
			double m[2];
			f[0] = EzIntFunction(0, it, r, x);
			m[0] = MsArray[0];

			for(int i=1; i<MsCount; i++)
			{
				//fill arrays
				f[1] = EzIntFunction(i, it, r, x);
				m[1] = MsArray[i];

				s += 0.5*(f[0]+f[1])*(m[1]-m[0]);	

				//update f[0]
				f[0] = f[1];
				m[0] = m[1];			
			}			
		}

		//push solvation
		s *= step;
		s /= -2*PI;
		ez.push_back(s);
	}

	
	//EzCache.PushV(iX, iY, iZ, ez);
	TimeInterpolation(ez, oResV, iTimes, iTimeCount);
	return true;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool LineAbDipole::DBzDt(double iX, double iY, double iZ, double *oResV, double *iTimes, int iTimeCount)
{
	if(!Layers) return false;
	if(!V || NeedRecalc) CreateX(iZ);

	std::vector<double> bz;
	bz.reserve(TimesCount);

	double r = sqrt(iX*iX + iY*iY);


	int n = GetSegmentsNum(iX, iY);
	double step = Length/n;
	BesselFunctions& bf = BesselFunctions::GetInstance();

	for(int it = 0; it<TimesCount; it++)
	{
		double s = 0;
		for(int sn = 0; sn<n; sn++)
		{
			double sx = -Length/2 + step * sn + step/2;
			double x = iX + sx;
			double r = FastMath::Sqrt(x*x + iY*iY);


			//numerical			
			double f[2];
			double m[2];
			f[0] = DBzDtIntFunction(0, it, r, x, iY);
			m[0] = MsArray[0];

			for(int i=1; i<MsCount; i++)
			{
				//fill arrays
				f[1] = DBzDtIntFunction(i, it, r, x, iY);
				m[1] = MsArray[i];

				BesselFunctions& b = BesselFunctions::GetInstance();
				s += step*(f[0]*b.J1(r*m[0]) + f[1]*b.J1(r*m[1]))*0.5*(m[1]-m[0]);//J1Int(f[0], f[1], r, m[0], m[1]);

				//update f[0]
				f[0] = f[1];
				m[0] = m[1];			
			}			
		}

		
		//push solvation
		s *= MU0/(2*PI);
		bz.push_back(s);
	}


	TimeInterpolation(bz, oResV, iTimes, iTimeCount);

	return true;
}


//--------------------------------------------------------------------------------------------------------------------------------------
bool LineAbDipole::By(double iX, double iY, double iZ, double *oResV, double *iTimes, int iTimeCount)
{
	if(!Layers) return false;
	if(!V || NeedRecalc) CreateX(iZ);

	std::vector<double> ez;
	ez.reserve(TimesCount);


	//try to extract from cache
	// 	if(EzCache.GetV(iX, iY, iZ, 0.01*r, ez))
	// 	{
	// 		TimeInterpolation(ez, oResV, iTimes, iTimeCount);
	// 		return true;
	// 	}


	int n = GetSegmentsNum(iX, iY);
	double step = Length/n;

	for(int it = 0; it<TimesCount; it++)
	{
		double s = 0;
		for(int sn = 0; sn<n; sn++)
		{
			double sx = -Length/2 + step * sn + step/2;
			double x = iX + sx;
			double r = FastMath::Sqrt(x*x + iY*iY);


			//numerical			
			double f[2];
			double m[2];
			f[0] = ByIntFunction(0, Layers[0].Sigma, it, r, x, iY);
			m[0] = MsArray[0];

			for(int i=1; i<MsCount; i++)
			{
				//fill arrays
				f[1] = ByIntFunction(i, Layers[0].Sigma, it, r, x, iY);
				m[1] = MsArray[i];

				s += 0.5*(f[0]+f[1])*(m[1]-m[0]);	

				//update f[0]
				f[0] = f[1];
				m[0] = m[1];			
			}			
		}

		//push solvation
		s *= step;
		s /= -MU0/(2*PI);
		ez.push_back(s);
	}


	//EzCache.PushV(iX, iY, iZ, ez);
	TimeInterpolation(ez, oResV, iTimes, iTimeCount);
	return true;
}


//--------------------------------------------------------------------------------------------------------------------------------------
LineAbDipole::~LineAbDipole()
{
	if(TimesArray) delete [] TimesArray;
}
